module polyroot
    implicit none
    integer,parameter :: N=3
    real :: poly(N+1)=(/-6.0,11.0,-6.0,1.0/),step=1d-5,error=1d5,rang=0.5
    real,save :: limit,sta
contains
    function f(p,x)
        implicit none
        real :: p(:),f,x
        integer :: i,M
        M=size(p)
        f=0.0
        do i=1,M
            f=f*x+poly(M-i+1)
        end do
    end function
    subroutine dp(p,q)
        real :: p(:),q(:)
        integer :: i
        q=p(2:)
        forall(i=1:size(q))q(i)=q(i)*i
    end subroutine
    real function ddvalue(p,q,r,x)
        real :: p(:),q(:),r(:),x
        real :: pv,qv,rv
        pv=f(p,x)
        qv=f(q,x)
        rv=f(r,x)
        ddvalue=(rv*pv-qv*qv)/pv**2.0
    end function
    subroutine root(r)
        implicit none
        real :: r(N),dd
        real :: ab(N),an,dpoly(N),ddpoly(N-1)
        integer :: i
        call dp(poly,dpoly)
        call dp(dpoly,ddpoly)
        ab=abs(poly(1:N))
        an=poly(N+1)
        limit=sum(ab)/abs(an)
        sta=-limit
        do i=1,N
            do while(sta<=limit)
                dd=ddvalue(poly,dpoly,ddpoly,sta)
                if(dd>error)then
                    exit
                end if
                sta=sta+step
            end do
            r(i)=sta
            sta=sta+rang
        end do
    end subroutine
end module
program main
    use polyroot
    implicit none
    real :: r(N)
    call root(r)
    write(*,*)r
end program
